/**
 * Code that implements the methods of the state class
 *
 * @author Carlos Colmenares (06-39380)
 * @author Kelwin Fernandez (07-40888)
 */

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <iostream>
#include "utils.h"
#include "state.h"

using namespace std;

int state::criterion_id;

state::state(clustering_instance* p, int nclusters){
    
    num_clusters = nclusters;
    prob = p;
    
    centroids = (double**)malloc(sizeof(double*)*num_clusters);
    closest_centroid = (int*)malloc( sizeof(int)*prob->n_data );
    distance_to_closest_centroid = (double*)malloc(sizeof(double)*prob->n_data);
    
    for (int c=0; c < nclusters; c++ ){
        centroids[c] = (double*)malloc(
                sizeof(double)*p->get_num_attributes());
    }
    
    initial_state();
    
    value = calc_criterion();
}

state::state(state* other){
    
    num_clusters = other->num_clusters;
    prob = other->prob;
    
    centroids = (double**)malloc(sizeof(double*)*num_clusters);
    closest_centroid = (int*)malloc( sizeof(int)*prob->n_data );
    distance_to_closest_centroid = (double*)malloc(sizeof(double)*prob->n_data);
    
    for (int c=0; c < num_clusters; c++ ){
        
        centroids[c] = (double*)malloc(
                sizeof(double)*prob->n_attributes);
        
        for (int a=0; a < prob->n_attributes; a++){
            centroids[c][a] = other->centroids[c][a];
        }
    }
    
    value = calc_criterion();
}

state::state(clustering_instance* p, vector<int> asign[]){
    
    int nclusters = 0;
    for (int i=0; i<p->n_data; i++){
        if ( asign[i].size() > 0 ){
            nclusters++;
        }
    }
    
    num_clusters = nclusters;
    prob = p;
    
    centroids = (double**)malloc(sizeof(double*)*num_clusters);
    closest_centroid = (int*)malloc( sizeof(int)*prob->n_data );
    distance_to_closest_centroid = (double*)malloc(sizeof(double)*prob->n_data);
    
    for (int c=0; c < nclusters; c++ ){
        centroids[c] = (double*)malloc(
                sizeof(double)*p->get_num_attributes());
    }
    
    nclusters = 0;
    for (int i=0; i<p->n_data; i++){
        if ( asign[i].size() > 0 ){
            
            for (int at=0; at<p->n_attributes; at++){
                centroids[nclusters][at] = 0;
            }
            
            for(int d=0; (unsigned int)d < asign[i].size(); d++){
                for (int at=0; at<p->n_attributes; at++){
                    centroids[nclusters][at] += p->data[ asign[i][d] ][at];
                }
            }
            
            for (int at=0; at<p->n_attributes; at++){
                centroids[nclusters][at] /= asign[i].size();
            }
            
            nclusters++;
        }
    }
    
    value = calc_criterion();
}

void state::filter_unused_clusters(){
    
    int* num_data_cluster = (int*)malloc( num_clusters*sizeof(int));
    memset(num_data_cluster, 0, num_clusters*sizeof(int));
    
    // Accumulate the mean distance of the elements in each centroid to itself
    for (int d=0; d < prob->n_data; d++){
        double min_dist = -1.0;
        int min_cluster = -1;
        
        // Compute the distance of each observation to its centroid.
        for (int c=0; c<num_clusters; c++){
            double dist = 
              euclidean_distance(prob->n_attributes, prob->data[d],
                                 centroids[c]);
            
            if ( min_cluster == -1 || dist < min_dist ){
                min_cluster = c;
                min_dist = dist;
            }
            
            closest_centroid[d] = min_cluster;
            distance_to_closest_centroid[d] = min_dist;
        }
        
        num_data_cluster[ min_cluster ]++;
    }
    
    int new_n_clusters = 0;
    for (int c=0; c<num_clusters; c++){
        if ( num_data_cluster[c] > 0 ){
            new_n_clusters++;
        }
    }
    
    double** centr = (double**)malloc(sizeof(double*)*new_n_clusters);
    int nc = 0;
    for(int c=0; c<num_clusters; c++){
        if ( num_data_cluster[c] != 0 ){
            centr[nc] = centroids[nc];
            nc++;
        }
    }
    
    free( centroids );
    centroids = centr;
    num_clusters = new_n_clusters;
    
    free(num_data_cluster);
    
    value = calc_criterion();
}

void state::initial_state()
{
    int next_point;
    
    //srand ( time(NULL) );
    srand ( rand() );
    for (int c=0; c < num_clusters; c++){
        next_point = rand() % prob->n_data;
        for (int attrib = 0; attrib < prob->n_attributes; attrib++){
            centroids[c][attrib] = prob->data[next_point][attrib];
        }
    }
}

state::~state() {
    for (int c=0; c < num_clusters; c++){
        free( centroids[c] );
    }
    free(centroids);
    free( distance_to_closest_centroid );
    free( closest_centroid );
}

state* state::get_neighbor(int cent, int point)
{
    state *neighbor = new state(this);
    
    for (int at = 0; at < prob->n_attributes; at++){
        neighbor->centroids[cent][at] += prob->data[point][at];
        neighbor->centroids[cent][at] /= 2.0;
    }
    
    neighbor->value = neighbor->act_criterion(cent);
    return neighbor;
}

state* state::get_random_neighbor(){
    double rnd = rand() / (double(RAND_MAX)+1);
    int cent = int(num_clusters * rnd);

    rnd = rand() / (double(RAND_MAX)+1);
    int point = int(prob->n_data * rnd);

    return get_neighbor(cent, point);
}

state* state::get_smooth_neighbor(int cent, int point)
{
    state *neighbor;
    
    int num_elements = 0;
    
    for (int d=0; d<prob->n_data; d++){
        if ( closest_centroid[d] == cent ){
            num_elements++;
        }
    }
    
    if ( num_elements != 0 ){
        
        neighbor = new state(this);
        
        for (int a=0; a<prob->n_attributes; a++){
                neighbor->centroids[cent][a] = 0.0;
        }
        
        for (int d=0; d<prob->n_data; d++){
            if ( closest_centroid[d] == cent ){
                for (int a=0; a<prob->n_attributes; a++){
                    neighbor->centroids[ closest_centroid[d] ][a] += prob->data[d][a];
                }
            }
        }
        for (int a=0; a<prob->n_attributes; a++){
            neighbor->centroids[cent][a] /= num_elements;
            neighbor->centroids[cent][a] += prob->data[point][a];
            neighbor->centroids[cent][a] /= 2;
        }
    } else {
        neighbor = get_neighbor(cent, point);
    }
    
    neighbor->value = neighbor->act_criterion(cent);
    return neighbor;
}

state* state::disturb(){
    
    state *sol = new state(this);
    int cent, point;
    
    srand ( time(NULL) );
    cent = rand() % num_clusters;
    
    int n_perturb = (num_clusters + 1)/2;
    while( n_perturb--){
        for (int a = 0; a < prob->n_attributes; a++){
            srand ( time(NULL) );
            point = rand() % prob->n_data;
            sol->centroids[cent][a] = prob->data[point][a];
        }
    }
    
    return sol;
}

double state::get_value(){
    return value;
}

bool state::acceptance_criteria(state* current, state* neighbor){
   
    switch(criterion_id){
        case 0:
            return current->value > neighbor->value;
        case 1:
            return current->value < neighbor->value;

        default:    cerr << "Wrong Criterion Asked!";
                    return false;

    
    }
}

double state::davies_bouldin_index(){
    double db = 0.0;
    
    // TODO: Este espacio de memoria se podria reservar globalmente para
    //       no estar pidiendolo a cada rato, revisar!
    double* sigma = (double*)malloc( num_clusters*sizeof(double));
    int* num_data_cluster = (int*)malloc( num_clusters*sizeof(int));
    
    memset(sigma, 0, num_clusters*sizeof(double));
    memset(num_data_cluster, 0, num_clusters*sizeof(int));
    
    // Accumulate the mean distance of the elements in each centroid to itself
    for (int d=0; d < prob->n_data; d++){
        double min_dist = -1.0;
        int min_cluster = -1;
        
        // Compute the distance of each observation to its centroid.
        for (int c=0; c<num_clusters; c++){
            double dist = 
              euclidean_distance(prob->n_attributes, prob->data[d],
                                 centroids[c]);
            
            if ( min_cluster == -1 || dist < min_dist ){
                min_cluster = c;
                min_dist = dist;
            }
            
            closest_centroid[d] = min_cluster;
            distance_to_closest_centroid[d] = min_dist;
        }
        
        sigma[ min_cluster ] += min_dist;
        num_data_cluster[ min_cluster ]++;
    }
    
    // If the centroid has elements, compute the mean distance.
    for (int c=0; c<num_clusters; c++){
        if (num_data_cluster[c] != 0){
            sigma[c] /= num_data_cluster[c];
        }else{
            sigma[c] = 0.0;
        }
    }
    
    // Calculate the Davies-Bouldin index
    
    for (int i=0; i<num_clusters; i++){
        
        // This inicialization is valid if the distance inter-cluster is
        // symmetric.
        for (int j=i+1; j<num_clusters; j++){
            
            double inter_cluster_distance =
                euclidean_distance(prob->n_attributes,
                                   centroids[i], centroids[j]
                                  );
            
            if ( inter_cluster_distance != 0 ){
                db = max(db, (sigma[i] + sigma[j])/inter_cluster_distance);
            }
        }
    }
    
    db /= num_clusters;
    
    free(sigma);
    free(num_data_cluster);
    
    return db;
}

double state::calc_criterion(){
    switch( state::criterion_id ){
        case 0: return davies_bouldin_index();
        case 1: return dunn_index();
        default: cerr << "Asked for an unexisting criterion";
                 return -1;
    }
}

double state::act_criterion(int centroid_moved){
    switch( state::criterion_id ){
        case 0: return davies_bouldin_index(centroid_moved);
        case 1: return dunn_index(centroid_moved);
        default: cerr << "Asked for an unexisting criterion";
                 return -1;
    }
}

double state::davies_bouldin_index(int centroid_moved){
    double db = 0.0;
    
    // TODO: Este espacio de memoria se podria reservar globalmente para
    //       no estar pidiendolo a cada rato, revisar!
    double* sigma = (double*)malloc( num_clusters*sizeof(double));
    int* num_data_cluster = (int*)malloc( num_clusters*sizeof(int));
    
    memset(sigma, 0, num_clusters*sizeof(double));
    memset(num_data_cluster, 0, num_clusters*sizeof(int));
    
    // Accumulate the mean distance of the elements in each centroid to itself
    for (int d=0; d < prob->n_data; d++){
        
        double min_dist = -1.0;
        int min_cluster = -1;
            
        if ( centroid_moved == closest_centroid[d] ){
            
            double dist = euclidean_distance(prob->n_attributes, prob->data[d],
                                        centroids[centroid_moved]);
            
            if ( dist <= distance_to_closest_centroid[d] ){
                distance_to_closest_centroid[d] = dist;
            } else {
                
                // Compute the distance of each observation to its centroid.
                for (int c=0; c<num_clusters; c++){
                    
                    if ( c == closest_centroid[d] )
                        continue;
                    
                    dist = 
                        euclidean_distance(prob->n_attributes, prob->data[d],
                                            centroids[c]);
                
                    if ( min_cluster == -1 || dist < min_dist ){
                        min_cluster = c;
                        min_dist = dist;
                    }
                }
                
                closest_centroid[d] = min_cluster;
                distance_to_closest_centroid[d] = min_dist;
            }
        } else {
            
            double dist = euclidean_distance(prob->n_attributes, prob->data[d],
                                        centroids[centroid_moved]);
            
            if ( dist < distance_to_closest_centroid[d] ){
                closest_centroid[d] = centroid_moved;
                distance_to_closest_centroid[d] = dist;
            }
        }
        
        sigma[ closest_centroid[d] ] += distance_to_closest_centroid[d];
        num_data_cluster[ closest_centroid[d] ]++;
    }
    
    // If the centroid has elements, compute the mean distance.
    for (int c=0; c<num_clusters; c++){
        if (num_data_cluster[c] != 0){
            sigma[c] /= num_data_cluster[c];
        }else{
            sigma[c] = 0.0;
        }
    }
    
    // Calculate the Davies-Bouldin index
    
    for (int i=0; i<num_clusters; i++){
        
        // This inicialization is valid if the distance inter-cluster is
        // symmetric.
        for (int j=i+1; j<num_clusters; j++){
            
            double inter_cluster_distance =
                euclidean_distance(prob->n_attributes,
                                   centroids[i], centroids[j]
                                  );
            
            if ( inter_cluster_distance != 0 ){
                db = max(db, (sigma[i] + sigma[j])/inter_cluster_distance);
            }
        }
    }
    
    db /= num_clusters;
    
    free(sigma);
    free(num_data_cluster);
    
    return db;
}

//This function is used by dunn_index
//It calculates the point p's closest centroid
//and actualizes the "closest centroid",
//distance_to_closest_centroid, max_intra_distance,
//and max_intra_d_point
void state::calc_closest_centroid(int p){

    closest_centroid[p] = 0;
    distance_to_closest_centroid[p] = 
        euclidean_distance(prob->n_attributes,
                prob->data[p], centroids[0]
                );

    for(int c=1;c<num_clusters;c++){
        double dist =
            euclidean_distance(prob->n_attributes,
                    prob->data[p], centroids[c]
                    );
        if( dist < distance_to_closest_centroid[p] ){
            distance_to_closest_centroid[p] = dist;
            closest_centroid[p] = c;
        }
    }

    //Always keep track of the maximum intra-distance:
    //The maximum distance between a point and its centroid.
    if( distance_to_closest_centroid[p] >
            max_intra_distance ){
        max_intra_d_point = p;
        max_intra_distance = distance_to_closest_centroid[p];
    }

}

//Among all clusters, calculate wich are the
//closest ones
void state::calc_closest_clusters(){

    closest_clusters = make_pair(0,1);
    min_inter_distance = 
        euclidean_distance(prob->n_attributes,
                centroids[0], centroids[1]
                );

    for(int i=0; i<num_clusters; i++){
        for(int j=i+1;j<num_clusters;j++){
            if(i==j)
                continue;

            double dist =
                euclidean_distance(prob->n_attributes,
                        centroids[i], centroids[j]
                        );
            if( dist < min_inter_distance ){
                min_inter_distance = dist;
                closest_clusters.first = i;
                closest_clusters.second = j;
            }
        }
    }
}

double state::dunn_index(){
    //For each point, calculate the centroid it should belong to
    max_intra_d_point = -1;
    max_intra_distance = 0;

    for(int p=0;p<prob->n_data;p++){
        calc_closest_centroid(p);        
    }

    //Calculate the minimum inter-distance: the two
    //closest centroids
    calc_closest_clusters();

    return min_inter_distance / max_intra_distance; 
}

double state::dunn_index(int centroid_moved){
    //If the centroid moved was the centroid with the
    //max intra distance, then everything will have
    //to be recalculated. You can optimize this.
    if( closest_centroid[max_intra_d_point] == centroid_moved ){
        return dunn_index();
    }

    //Check every point and recalculate the closest
    //centroid and it's distance
    for(int p=0;p<prob->n_data;p++){
        
        //If the moved centroid was the closest one to p,
        //recalculate completely the closest
        //centroid. Else, only calculate if the moved
        //centroid is closer
        if( closest_centroid[p] == centroid_moved ){
            calc_closest_centroid(p);
        }
        else {
            //New distance to centroid
            double dist = 
                euclidean_distance(prob->n_attributes,
                        prob->data[p], centroids[centroid_moved]
                        );
            if( dist < distance_to_closest_centroid[p] ){
                distance_to_closest_centroid[p] = dist;
                closest_centroid[p] = centroid_moved;

                //Always keep track of the maximum intra-distance:
                //The maximum distance between a point and its centroid.
                if( distance_to_closest_centroid[p] >
                        max_intra_distance ){
                    max_intra_d_point = p;
                    max_intra_distance = distance_to_closest_centroid[p];
                }

            }
        }
    }
    //Now, if the moved centroid was one of the closest
    //centroids, then recalculate everything, else
    //just if it's now closer to another point
    if( closest_clusters.first == centroid_moved ||
        closest_clusters.second == centroid_moved ){
        
        calc_closest_clusters();

    } else {
        int i = centroid_moved;
        for(int j=0;j<num_clusters;j++){
            if(i==j)
                continue;

            double dist =
                euclidean_distance(prob->n_attributes,
                        centroids[i], centroids[j]
                        );
            if( dist < min_inter_distance ){
                min_inter_distance = dist;
                closest_clusters.first = i;
                closest_clusters.second = j;
            }
        }
    }

    return min_inter_distance / max_intra_distance; 
}

void state::assign_nearest(){
    
    for (int d=0; d < prob->n_data; d++){
        calc_closest_centroid(d);
    }

}

void state::recompute_centroids(){
    
    int* num_elements = (int*)malloc(num_clusters*sizeof(int));
    memset( num_elements, 0, num_clusters*sizeof(int) );
    
    for (int d=0; d<prob->n_data; d++){
        num_elements[ closest_centroid[d] ]++;
    }
    
    for (int c=0; c<num_clusters; c++){
        if ( num_elements[c] != 0 ){
            for (int a=0; a<prob->n_attributes; a++){
                centroids[c][a] = 0.0;
            }
        }
    }
    
    for (int d=0; d<prob->n_data; d++){
        for (int a=0; a<prob->n_attributes; a++){
            centroids[ closest_centroid[d] ][a] += prob->data[d][a];
        }
    }
    
    for (int c=0; c<num_clusters; c++){
        if ( num_elements[c] != 0 ){
            for (int a=0; a<prob->n_attributes; a++){
                centroids[c][a] /= num_elements[c];
            }
        }
    }
    
    calc_criterion();
    
    free( num_elements );
}

void state::print(){
    cout << value << endl;
};
